Non-quantum liquefaction of coherent gases 
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We show that a gas-to-Uquid phase transition at zero temperature may occur in a coherent gas 
of bosons in the presence of competing nonlinear effects. This situation can take place both in 
atomic systems like Bose-Einstein Condensates in alkalii gases with two and three-body interactions 
of opposite signs, as well as in laser beams which propagate through optical media with Kerr 
(focusing) and higher order (defocusing) nonlinear responses. The liquefaction process takes place 
in absence of any quantum effect and can be formulated in the frame of a mean field theory, in 
terms of the minimization of a thermodynamic potential. We also show numerically that the effect 
of linear gain and three-body recombination also provides a rich dynamics with the emergence of 
self-organization behaviour. 
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I. INTRODUCTION 

Paraxial propagation of linearly polarized laser beams 
through transparent optical media with intensity depen- 
dent refractive index is mathematically equivalent to the 
free evolution of the wavefunction order parameter used 
in the mean-field description of a two-dimensional gas 
of N interacting atoms in a Bose-Einstein Condensate 
(BEC) at temperature T = OK Both systems can be 
modeled by identical nonlinear Schrodinger equations [2]. 
For photons in the laser beam, the x*^") component of the 
nonlinear optical susceptibility plays the same role as n- 
body interactions between atoms in the cloud and the 
propagation constant can be identified with a chemical 
potential for the light distribution. As all the photons in 
a coherent wave are equal, the laser beam can be treated 
on equal foot as any system of N identical interacting 
bosons at zero temperature |3j. 

The previous point of view, which takes into account 
the equivalence between laser beams and BECs of ultra- 
cold atoms, has led to an interesting suggestion made by 
Chiaojl], who recently proposed to verify the superfluid- 
ity of coherent light, in analogy with degenerate quan- 
tum atomic gases. More recently, similar concepts have 
been successfully used to analyze condensation phenom- 
ena of nonlinear waves fSI and quantum phase transitions 
of photons in periodic lattices [B]. 

In Chiao's model the key point of the analysis is "same 
equations, same predictions" and therefore photons from 
a monochromatic laser source are considered as an ideal 
bosonic gas at zero temperature in which continuous 
phase transitions (CPT) can take place due to long range 
quantum fluctuations around the ground state [7]. These 
critical phenomena are thus called quantum phase tran- 
sitions (QPT)|j8 to distinguish them from the standard 
phase changes which are well-known in classical thermo- 
dynamics. 

As we will show in this work, CPT may occur in any 
classical system at zero temperature without long range 
quantum correlations involved, if opposite nonlinear in- 



teractions are present. In the case of Chiao's "superfluid 
light" the phase transition is produced by effect of a defo- 
cusing intensity-dependent refractive index [5] and thus a 
waveguide is used to avoid spreading of the beam (in the 
same way as magneto-optical traps are employed to hold 
atomic BECs). However, as superfluidity may occur both 
in gases and liquids, it cannot be considered as a trace of 
the presence of a liquid state[Tni E]; it is also required 
the appearance of surface tension effects[T2]. Moreover, 
in Chiao's model, the interactions between particles are 
repulsive and they cannot drive a gas-liquid transition. 

Thus, in this work we will follow the same lines of 
thought to suggest the possibility of obtaining a gas- 
to-liquid phase transition in a classical gas at T = OK 
described in a mean field theory by the so-called cubic- 
quintic (CQ) model with competing nonlinearities. Sev- 
eral pioneering works have highlighted the interesting 
properties of this CQ-model|13J. Cavitation, superflu- 
idity and coalescence have been investigated [13] TF in 
the context of liquid He, where the model is a simple 
approach if nonlocal interactions are not taken into ac- 
count. Stable optical vortex solitons and the existence of 
top-flat states have been also reported in optical mate- 
rials with CQ optical susceptibility [THJ [T7]. The surface 
tension properties that appear in this system[T8] have 
been considered as a trace of a "liquid state of light" [TD] . 

On the other hand, recent experiments about fllamen- 
tation of high-power laser pulses in CS2 have shown that 
the CQ nonlinearity is achievable in this material[19 . It 
has been also suggested that atomic coherence may be 
used to induce a giant CQ-like refractive index of a Rb 
gasfn . Thus, the practical realization of the flrst "liq- 
uid of light" state as an example of non-quantum lique- 
faction at zero temperature is close. In BEC systems, a 
CQ-model can be used in the mean field description of 
an ultracold gas at zero temperature in the presence of 
Efimov states with tunable two- and three-body interac- 
tions, which have been recently proposed |20j. 
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II. MATHEMATICAL ANALYSIS 

The cubic-quintic model describes a coherent bosonic 
system of N particles with two- and three-body interac- 
tions. The mathematical formulation of the mean field 
theory yields a generalised non- linear Schrodinger (NLS, 
also called Gross-Pitaevskii) equation of the form: 



.5* 



-Vl'f + 7|^'|2^' - S\^\^^ = 0, 



(1) 



If the system modeled by the previous equation is a 
photon gas, the above NLS describes paraxial prop- 
agation of a continuous linearly-polarized laser beam 
of wavelength A in a nonlinear medium with a refrac- 
tive index depending on the intensity / in the form 
n — uq + — n^I^ and the adimensional variables 
are: rj the propagation distance multiplied by 27r/A, j^fp 
the beam irradiance multiplied by the Kerr coefficient 
712 = 7, 714 = 5 an adimensional parameter indicating 
the strength of the quintic nonlinear optical susceptibil- 
ity, and V^^ = / dx^ + / dy^ the transverse Laplacian 
operator, where x and y are the transverse spatial dimen- 
sions multiplied by 27rY^2no/A. 

In the case of a two-dimensional atomic BEG tightly 
trapped along one axis by a parabolic potential of fre- 
quency v±^ and thickness = ^hjrav^, the adimen- 
sional variables correspond to: r\ the time in units 
of v'^ , l^fp the atomic density multiplied by two- 
body coefficient 7, (5 an adimensional parameter indi- 
cating the strength of the three-body interactions, and 
— /dx^ + jdy^ the transverse Laplacian opera- 
tor, where x and y are the transverse spatial dimensions 
divided by r^. 

Recent experiments for beam propagation in C&llQj 
at A = 800nm yield the following values for the above 
parameters in the case of laser beams: — 1.6, ?t-2 = 
3 • l^-^^cm^jW and n4 = 2 • IQ-^^cm^/VF^ Other ma- 
terials like air [21 or chalcogenide glasses [22] seem to dis- 
play the C-Q behavior usually accompanied by ionization 
and non-linear losses. It has been also pointed the pos- 
sibility of engineering this type of optical response by 
quantum techniques which allow to access this nonlin- 
ear regime with miliwatt ultrastabilized lasers [11]. For 
atomic BEG systems it has been proposed that a com- 
bination of two-body (attractive) and three-body (re- 
pulsive) elastic interactions can yield liquid hoselets^I^. 
This behaviour can be explained in terms of the Efimov 
states [2(1. However, three-body scattering in BEGs has 
inelastic contributions and yields highly nonlinear losses. 
This means that in the most general case the coefficient 
(5 in equation ([T]) may be complex both for laser systems 
as well as for atomic gases. 

The above GQ-NLSE admits soliton-like solutions of 
finite size[T3] of the form ^'^(a;, y, 77) — A{x,y)e~'^'^^ , 
being /i the propagation constant in the case of light 
and the chemical potential for atomic BEGs. These 
solitons can only be calculated numerically and coex- 
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FIG. 1: [Color online] Maximum squared amplitude of differ- 
ent types of stationary solutions of Eq. ([l|. The continuous 
and dotted curves correspond to plane waves. The dashed 
line represents numerically calculated localized eigenstates. 
The continuous and dashed curves join at 7|^maa:|^ ~ 0.75, 



Horizontal lines indicates the critical values 



7j^mai:|^ = 0.5 and 7|Amaa;P ~ 0.75. Labeled points cor- 
respond to solitons that are considered as special examples in 
the text. 



ist with plane waves solutions of constant amplitude 
A{x,y,r]) = Ae"***'', which lead by substitution in 
Eq.Q to /i = 5\A\'^ ~ 7|A|2. In Fig (l] we have plotted 
the maximum value of 7|5'P vs. /i in units of jieq = 
—0.18757^/(5, for different kinds of stationary solutions 
of Eq. ([l]) . The continuous and dotted lines correspond 
respectively to stable and modulationally unstable plane 
waves [29], whereas the dashed line stands for numerically 
calculated localized eigenstates [2S]. It is known [TO that 
the shape of the solitons (see dashed profiles in Figj2]) 
vary from quasi-gaussian shapes for low powers to almost 
square profiles for beam amplitudes close to a certain 
critical value of fj, — ^eq [ISl Hi] • At this point the size of 
the solutions tends to infinity whereas the amplitude of 
the beam stabilizes at \Amax\^ = \Aeq\^ = 0.757/(5. We 
will now clarify the reasons of this behavior by using a 
thermodynamic model. 



III. 



THERMODYNAMIC MODEL 



The appropriate tool to study the equilibrium con- 
dition as a function of the number of particles N = 
/ Y^\^dxdy is Landau's grand potential 



n = H - nN = 



(2) 



dxdy 



1*1 



1*1 



*l 



where the Lagrange multiplier /i is the chemical poten- 
tial for a BEG or the propagation constant in the case 
of an optical system. For our two-dimensional model the 
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FIG. 2: [Color online] Pressure plots of the plane wave so- 
lution branches of Fig[l] (green solid, stable; purple dashed, 
unstable) and pressure at the center [x — 0) of the eigen- 
states (thick dotted line) corresponding to different points in 
the dashed line of Fig [l] Insets; Profiles of the eigenstates 
(dashed) and pressure distribution (solid). The x-axis indi- 
cates the transverse size of the numerically calculated pressure 
distributions of different eigenstates of eq.([T]). We have mul- 
tiplied both 7|^P and p by 10 and 100 respectively in a) and 
p by 10 elsewhere. 



area is S and we have that p = —dQ,/dS is the pres- 
sure. The equilibrium configurations are obtained when 
U, is minimized. In particular dfl/dS = implies that 
the pressure has to be zero, which for the plane waves 
of Fig[l| yields dfl/dS = -^\A\^ + f - fi\A\^ = 0. 
There are two possibilities: the trivial case \Aeq\ = and 

a uniform phase, with /ieq — — fgy and l^egP — 0.757/(5. 
In the language of field theory, these solutions are the 
two possible "vacuum" states of the system. As it can 
be seen in Fig. [l] the two vacua can be "connected" by 
the soliton solutions of the dashed line in Fig. [l] which 
constitute the instantons of the theory It is interest- 
ing to notice that the non-zero vacuum solution implies 
a spontaneously symmetry breaking of the global phase 
symmetry, which is preserved in the case with \A\ = 0. 

This result about the pressure has already been dis- 
cussed by authors in ref [14| within the framework of the 
Madelung Transformation (MT), which allows one to es- 
tablish a formal analogy between nonlinear optics and 
classical fluid dynamics [27] • With the aid of the MT, 
an analytical expression for the effective pressure has 
been derived [H] . Very remarkably, this effective pres- 
sure, which is an approximation since it does not take 
into account the quantum-mechanical pressure term[27], 
calculated for the plane-wave solutions of the CQ-NLSE 
has the same mathematical expression as the one shown 
before. Unlike the previous works related in the refer- 
ences, we have calculated the effective pressure by means 
of the potential given by eq.([2]) which contains all the 
relevant information about the nonlinear system which 



is currently being studied. 

While the analysis of the inflnite plane-wave solutions 
has already been considered in the literature |14j , our dis- 
cussion of the pressure also applies to the localised soli- 
ton solutions, being a central contribution of the present 
work. We also flnd a signiflcant difference between the 
low-power quasi-gaussian solutions and the high-power 
top-flat solutions of Eq.([T]). The expression for plane- 
waves and top-flat eigenstates can be approximated by 
the same "reduced" expression for p given by both the 
MT and the omega density calculation. In fact, we have 
shown in Figjl] that the two branches of solutions, i.e. 
the localised solitons and the stable plane waves, merge 
as fi approaches ^eq, having a radius increasing to infin- 
ity, thus in this case it is justified to neglect |V_l^'P in 
Eq.(|2]). On the other hand, for the quasi-gaussian eigen- 
states the previous gradient term is not negligible. In 
this case, we should calculate p numerically in the lack 
of analytical solutions. 

The analysis of the pressure can be refined calculat- 
ing numerically its distribution for different eigenstates 
as it is plotted in Figj2] The pressure of the plane wave 
branches in Figjl] is plotted as a function of /i//ieq. Very 
remarkable is the fact that the pressure of the stable 
branch is higher than its counterpart of the unstable 
branch. This implies that the stable plane waves free 
energy density is smaller than the free energy of the mod- 
ulationally unstable solutions branch. In Figj2]it is also 
plotted the curve of the eigenstates central pressure. As 
it can be seen in the graph, the curve corresponding to 
the eigenstates is bounded by the two curves of the plane- 
waves branches, so that the existence domain of the flla- 
ments phase is limited by them. 

In the insets of Figj2] we show the shape profiles of 
the eigenstates (dashed line) superposed with their effec- 
tive pressure profiles (solid lines) which have been conve- 
niently rescaled to fit in the graph. As it can be appre- 
ciated in inset a), solitons with <C |A*eq| have smooth 
pressure distributions with a central maximum located 
at the centroid of the soliton and two negative-valued 
minima. As the value of increases, the soliton profiles 
and their corresponding shapes of p narrow, reaching a 
minimum width at fi/fieq — 0.5, which corresponds to 
a filament soliton solution with the same peak ampli- 
tude as the "critical" plane wave with \A\ = \Acr\- This 
plane-wave marks the border between stable and modu- 
lationally unstable plane waves When the chemical 
potential reaches the value |^| — \fifu\ (see inset b), the 
absolute maximum of p for eigenstates is obtained, i.e, 
this filament has the minimum of dfl/dS at its centroid. 
We also consider that close to = |/iHg|, the liquid top- 
flat eigenstates begin to exist [TO . In fact, the curve of 
stable localised solutions in Fig|T]has a maximum in that 
region and both localised solution and stable plane- waves 
branches seem to merge there. For < < l/ieql, 

the pressure maximum is located in a flat region (see in- 
set c), tending to zero as approaches the critical value 
IfJieql, point d) in Fig|T| over which no localized solutions 
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exist HI]. 



IV. CONDENSATION IN THE PRESENCE OF 
LINEAR GAIN AND NONLINEAR LOSSES 

In this section, we will provide a set of numerical simu- 
lations showing the condensation process, i.e., the phase 
transition from the gaseous phase to an homogeneus co- 
herent "liquid" plane wave solution corresponding to the 
upper branch in Figjl] In order to achieve this result, 
we will perturb an unstable plane wave, corresponding to 
the lower branch in Fig[T] with a randomly varying noise. 
This will produce a filamented phase, made of coherent 
structures whose shape remains qualitatively unchanged, 
up to smaller scale fluctuations [231. 

In order to achieve the non-quantum liquefaction we 
have included both a linear incoherent pumping mecha- 
nism and nonlinear three-body losses [ST]. In other words, 
we have considered 6 — S^+iS^ and we have introduced a 
linear gain term iT'i' in eq. ([l]) . This corresponds to a con- 
tinuous load of particles in the system. Note that in this 
way, we are describing a more realistic non-conservative 
version of eq.Q, which models an experimentally achiev- 
able scenario in the framework of current BEC experi- 
ments. In fact, in condensed matter systems it is possible 
to control the load of particles and two and three-body 
recombination within the coherent atomic cloud 1 20 



Ginzburg-Landau (G-L) equations used to describe wide- 
aperture laser cavities [33] . Although it is always possible 
to control the linear gain introduced in the system, three- 
body losses are often imposed by the nonlinear response 
of the material, so it is not possible to manage the dissipa- 
tion terms of the system. However, by means of electro- 
magnetic induced transparency techniques, it is possible 
to customize the nonlinear optical response of cold atomic 
ensembles like Rb|TT] so that the nonlinear refractive in- 
dex corresponds to the one given by the modified eq.Q 
analyzed in the current section of the paper. Therefore 
our model, and its predicted phenomenon of liquefaction 
that we will demostrate below, can correspond both to 
realistic BEC and nonlinear optical systems. 

In Fig|3] the initial state consists of an incoherent set 
of filaments with a randomly varying phase distribution. 
Within this apparent disorder, some coherent uncorre- 
lated structures (filaments) exist and can be observed in 
Figj3]a). In this simulation, we have considered a non- 
linear parameter regime where the linear gain was not 
enough to compensate the three-body dissipative term. 
As a consequence, certain degree of coherence is lost since 
the filaments disappear and only the noisy background 
is observed, as shown by snapshot c) in Fig|3|. Starting 




FIG. 3: [Color online] Numerical simulation of the evolution 
of a set of filaments in presence of linear gain and three-body 
losses. Propagation algorithm parameteres: 7 = 1, F = 10~^, 
S = 1 + O.li. Top: pseudocolor maps of the amplitude. Bot- 
tom: pseudocolor maps of the phase corresponding to the 
amplitudes above. As it can be appreciated in the sequence 
of snapshots, the initial condition which has a certain level 
of organization, evolves towards a complete disordered situ- 
ation, i.e., the diluted disordered gaseous phase arises. The 
frames correspond to values of the adimensional variable rj: 
(a, d), 1000 (b, e) and 4000 (c, f). 

On the other hand, in nonlinear optics, this kind of 
nonlinear models are well-known in the frame of complex 



FIG. 4: [Color online] Same as in Figjs] but with a modified 
value of the linear gain F = 8 ■ 10^''. Other parameters as 
in Fig|3] In this simulation, the final state is a homogeneous 
background of coherent liquid with amplitude \A\ = [v4eq|, but 
it is remarkable that the wavefront presents some topological 
defects (vortices) which distort the phase distribution of the 
resulting liquid plane wave as it can be observed in the phase 
maps on the bottom row. The frames correspond to values of 
the adimensional variable rj: (a, d), 1000 (b, e) and 4000 (c, 
f). 

from the same initial condition but increasing the linear 
gain term over a certain threshold, we have performed 
the simulation shown in Fig|4] In this case, we see that 
the system evolves towards an homogeneus plane wave 
by the combined effect of adding particles to the initial 
random state and dissipation due to many-body inelastic 
processes. This relaxation process is well-known in the 
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context of the complex G-L equations and is attributted 
to the non-conservative nature of the model [32] • Never- 
theless, as our theory predicts, the system will tend to 
form a particular plane wave of the stable "liquid" upper 
branch of Figjl] the one with zero pressure. This illus- 
trates the tendency of the system to reach the non-zero 
vacuum state with \A\'^ ~ 3j/A6. Very remarkably, dur- 
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FIG. 5: [Color online] Numerical simulation of the evolution 
of a set of filaments in presence of C-Q nonlinear gain and 
loss. The linear gain term used in the previous simulations 
has been replaced here by a nonlinear gain term ixl^l^j where 
X = 0.1. Other parameters are: 7=1,5 = 1-1- O.li. Top: 
pseudocolor maps of the amplitude. Bottom: pseucolor maps 
of the phase corresponding to the amplitudes above. Quali- 
tatively, the same results of Fig|4] are obtained. However, it 
is important to note that in this case, the robust topologi- 
cal structures with nonzero vorticity which appear during the 
dynamical liquefaction process are annihilated, resulting in 
a final homogeneous background with amplitude \A\ — \Aeq\ 
and constant phase distribution, as it can be seen in snap- 
shot f). The frames correspond to values of the adimensional 
variable rj: (a, d), 1000 (b, e) and 4000 (c, f). 

ing the process pairs of vortex- antivortex with topolog- 
ical charges to„ = l,mav = — 1 are formed (see bottom 
phase maps in Figj4]), so that the constant phase of the 



emerging coherent wave remains hidden by the overlap- 
ping of the different vortex rotating phase-distributions. 
Vortices are very robust topological structures [T^] and 
in our simulations they remain stable as far as we could 
follow the numerical simulations. 

Finally, we have considered the effect of replacing the 
linear gain term by a nonlinear gain term of the form 
«x|^p. In this situation, the numerical results are quali- 
tatively the same as in Fig|4] although it can be observed 
that in the last stage of the evolution the vortices are an- 
nihilated. In snapshot f) of Figjl] it can be seen how 
the underlying plane wave, which emerges after the dy- 
namical process as described above, has a homogeneus 
phase distribution. We think that this result is a very 
interesting example of a self-organization process, where 
coherence is produced from disorder, with evident prac- 
tical applications. 

V. CONCLUSIONS 

We have shown in the present work that a system of 
N equal bosons in a system with competing nonlineari- 
ties can undergo a phase transition from a gas to a liquid 
state. The process takes place at zero temperature with- 
out any quantum effect and it is only ruled by nonlinear 
interactions. Finally, we have shown that a cubic-quintic 
medium with complex susceptibilities exhibiting linear 
gain and nonlinear losses, will tend to produce a homoge- 
neous phase liquid distribution starting from a collection 
of non-correlated filaments. This opens the door for ex- 
periments in the field of BEG systems in ultracold gases. 
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